*! Version 2.0.6  17may2006
* version 1 Marcelo Moreira and Brian Poi
* Estimated IV regression by 2sls or liml
* version 1.0.0 created 20021105 by Brian Poi
* version 2.0.0 created 20051001 by Anna Mikusheva
* version 2.0.1 20051205 code check and reformat for SJ inclusion by BPP
* version 2.0.2 changed final display format
* version 2.0.3 bug fixes
* version 2.0.4 bug fix regarding first-stage F statistic
* version 2.0.5 bug fix regarding multicollinearity
* version 2.0.6 renamed e(LR_p) -> e(p_LR) and likewise for AR, LM stats
program define condivreg, eclass byable(recall)

	version 8.2
	
	if replay() {
		if "`e(cmd)'" != "condivreg" {
			error 301
		}
		syntax [, level(integer -1) ]
		if `level' > -1 {
			di as error "cannot specify {cmd:level()} on replay"
			exit 198
		}
		myprint
		myprintregions
		exit
	}
		 
	/* First get the list of variables.	*/
	local n 0
	gettoken lhs 0 : 0, parse(" ,[") match(paren)
	IsStop `lhs'
	if `s(stop)' { 
		error 198 
	}
	while `s(stop)'==0 {
		if "`paren'"=="(" {
			local n = `n' + 1
			gettoken p lhs : lhs, parse(" =")
			while "`p'"!="=" {
				local end`n' `end`n'' `p'
				gettoken p lhs : lhs, parse(" =")
			}
			tsunab end`n' : `end`n''
			tsunab rinst : `lhs'
		}
		else {
			local exog `exog' `lhs'
		}
		gettoken lhs 0 : 0, parse(" ,[") match(paren)
		IsStop `lhs'
	}
	local 0 `"`lhs' `0'"'

	tsunab exog : `exog'
	tokenize `exog'
	local ry1 "`1'"
	local 1 " "
	local rexog `*'
	loc ry2 "`end1'"
	syntax [if] [in], [nocons noinstcons liml 2sls ar lm interval 	/*
		*/	   level(integer $S_level) test(real 0)]
	
	marksample touse
	markout `touse' `ry1' `ry2' `rexog' `rinst'

/* Some syntax checking.	*/
	if ("`liml'" != "" & "`2sls'" != "") {
		di as error "Only one of liml or 2sls is allowed."
		exit 198
	}
	if ("`liml'" != "") {
		local est_type "liml"
	}
	else {
		local est_type "2sls"
	}

	quietly {
/* check for multicollinearity */
		tempvar one  AAA  vv dd 
		loc k : word count `rinst'
		loc p : word count `rexog'
/* BPP added 2006-04-18 */
		_rmcoll `rexog' `rinst', `cons'
		local newcnt : word count `r(varlist)'
		local oldcnt : word count `rexog' `rinst'
		if `newcnt' != `oldcnt' {
			noi di as error "Multicollinearity!"
			exit 498
		}
/* End BPP addition */
		gen `one' = 1
	 	loc XZ "`rexog' `rinst'"
		if ("`cons'" == "") {
			mat accum `AAA' = `XZ'
			loc k = `k' + 1
		}
		else{
			mat accum `AAA' = `XZ', nocons
		}
/*
BPP commented out this block 2006-04-18
		mat symeigen `vv' `dd' = `AAA'
		if (`dd'[1,`k'+`p']<0.000000001) {
			noi di as error "Multicollinearity!"
			exit
		}
*/

		tempname bbb 
		sca `bbb'=`test'
		if ("`ar'"!="") { 
			loc AR "ar" 
		}
		if ("`lm'"!="") { 
			loc LM "lm" 
		}
 
/* Compute the 2sls and liml results.	*/
		tempname beta var rmse rss mss rsquared adjrsq fstat 	/*
			*/	xpx limlbeta
		myliml `ry1' `ry2' "`rexog'" "`rinst'" "`cons'" 	/*
			*/	"`instcons'" `touse'
		mat `beta' = r(beta)
		sca `limlbeta' = `beta'[1,1]

		if "`est_type'" == "2sls" {
			my2sls `ry1' `ry2' "`rexog'" "`rinst'" "`cons'"	/*
				*/	"`instcons'" `touse'  
		}
		mat `beta' = r(beta)
		mat `var' = r(var)
		if "`instcons'" == "" {
			mat colnames `beta' = `ry2' `rexog' _cons
			mat rownames `var' = `ry2' `rexog' _cons
			mat colnames `var' = `ry2' `rexog' _cons
								
		}
		else {
			mat colnames `beta' = `ry2' `rexog'
			mat rownames `var' = `ry2' `rexog'
			mat colnames `var' = `ry2' `rexog'
		}
		sca `rss' = r(rss)
/* Compute the general stuff that goes along. */
		tempvar junk
		gen double `junk' = `ry1'^2 if `touse'
		su `junk', meanonly
		loc capn = r(N)
		sca `mss' = r(sum) - `rss'
		if ("`cons'" == "") {
			su `ry1' if `touse', meanonly
			sca `mss' = `mss' - r(sum)^2 / `capn'
		}
		loc dfm : word count `rexog' `ry2'
		loc dfr = `capn' - `dfm'
		if ("`cons'" == "") {
			loc c = 1
			loc dfr = `dfr' - 1
		}
		else {
			loc c = 0
		}
		sca `rmse' = sqrt(`rss'/`dfr')
		sca `rsquared' = `mss' / (`mss' + `rss')
		sca `adjrsq' = 1 - (1-`rsquared')*(`capn'-`c') /	/*
			*/		(`capn'-`c'-`dfm')
		/* F statistic for 2SLS/Wald for LIML. */
		if "`cons'" != "" {
			sca `fstat' = .
		}
		else {
			tempname c dif vinv
			mat `c' = J(1, (`dfm' + 1), 0)
			su `ry1' if `touse', meanonly
			mat `c'[1, (`dfm'+1)] = r(mean) 
			mat `dif' = `beta' - `c'
			mat `vinv' = inv(`var')
			mat `fstat' = `dif'*`vinv'*`dif''
			sca `fstat' = trace(`fstat')
			if "`est_type'" == "2sls" {
				sca `fstat' = `fstat' / `dfm'
			}
		}

/* Get the first-stage stats.	*/
		tempname f1 r21 r2_a1 df_m1 df_r1
		if "`instcons'" == "" {
			reg `ry2' `rexog' `rinst'
		}
		else {
			reg `ry2' `rexog' `rinst', nocons
		}
		test `rinst'
		sca `f1' = r(F)
		sca `r21' = e(r2)
		sca `r2_a1' = e(r2_a)
		sca `df_m1' = r(df)
		sca `df_r1' = r(df_r)
		
/* Post results.  Use two different syntaxes of 2sls and liml
   so that we get t-stat for 2sls and z-stat for liml. */
   		tempvar touse2
   		gen byte `touse2' = `touse'	/* gets destroyed by eret post */
		if "`est_type'" == "2sls" {
			eret post `beta' `var', esample(`touse2') 	/*
				*/	depname("`ry1'") obs(`capn') dof(`dfr')
		}
		else {
			eret post `beta' `var', esample(`touse2') 	/*
				*/	depname("`ry1'") obs(`capn')
		} 
		eret local depvar "`ry1'"
		eret scalar df_m = `dfm'
		eret scalar df_r = `dfr'
		if "`est_type'" == "2sls" {
			eret scalar F = `fstat'
		}
		else {
			eret scalar wald = `fstat'
		}
		eret scalar H0_b = `bbb'
		eret scalar r2 = `rsquared'
		eret scalar r2_a = `adjrsq'
		eret scalar rmse = `rmse'
		eret scalar mss = `mss'
		eret scalar rss = `rss'
		eret scalar F_first = `f1'
		eret scalar df_m_first = `df_m1'
		eret scalar df_r_first = `df_r1'
		eret scalar r2_first = `r21'
		eret scalar r2_a_first = `r2_a1'
		eret scalar b_instd_liml = `limlbeta'
		eret scalar level = `level'
		eret local exog "`rexog'"
		eret local inst "`rinst'"
		eret local insts "`rexog' `rinst'"
		eret local instd "`ry2'"
		if "`est_type'" == "2sls" {
			eret local model "2sls"
		}
		else {
			eret local model "liml"
		}
		if "`instcons'" == "" {
			eret local instcons "yes"
		}
		else {
			eret local instcons "no"
		}
		if "`cons'" == "" {
			eret local cons "yes"
		}
		else {
			eret local cons "no"
		}

/* Common block for testing and conf sets */
 		if "`cons'" == "" {
			loc exog "`rexog' `one'"
		}
		else {
			loc exog "`rexog'"
		}
/* preserve currently posted results; -reg- will clobber them otherwise */
		tempname myest
		_estimates hold `myest'
		tempvar y1 y2
		if ("`exog'" != "") {
			foreach v in y1 y2 {
				reg `r`v'' `exog' if `touse', nocons
				predict double ``v'' if `touse', residuals
			}
		}
		else {
			gen double `y1' = `ry1' if `touse'
			gen double `y2' = `ry2' if `touse'
		}
/* Regress instruments on exog.	*/
		tempname ehold
		loc inst = ""
		loc n = 1
		foreach v in `rinst' {
			tempvar inst`n'
			if ("`exog'" != "") {
				reg `v' `exog' if `touse', nocons
				predict double `inst`n'' if `touse', residuals
			}
			else {
				gen double `inst`n'' = `v' if `touse'
			}
			loc inst "`inst' `inst`n''"
		}
/* Compute Omega. */
		tempname mzy1 mzy2 n k omega
		if "`instcons'" == "" {
			reg `y1' `inst' if `touse'
		}
		else {
			reg `y1' `inst' if `touse', nocons
		}
		predict `mzy1' if `touse', residuals
		if "`instcons'" == "" {
			reg `y2' `inst' if `touse'
		}
		else {
			reg `y2' `inst' if `touse', nocons
		}
		predict `mzy2' if `touse', residuals
		mat accum `omega' = `mzy1' `mzy2', noconstant
/* Safe to restore my results */
		_estimates unhold `myest'
		count if `touse'
		loc n = r(N)
		loc k : word count `inst'
		loc p : word count `exog'
		
		if ("`instcons'" == "") {
			loc k = `k' + 1
		}
		mat `omega' = `omega' / (`n'-`k'-`p')
		if (("`instcons'" == "")& ("`cons'" == "")) {
			loc df = `k' - 1	
		}
		else {
			loc df = `k' 
		}

/*------Confidence sets -------------*/

/*
Test		Result type		Interval
-----------------------------------------------------------------------
CLR		1			Empty set
		2			[x1, x2]
		3			(-infty, +infty)
		4		    (-infty, x1] U [x2, infty)
		
AR		1			Empty set
		2			[x1, x2]
		3			(-infty, +infty)
		4		    (-infty, x1] U [x2, infty)
		
LM		1			Not used (not possible)
                2			[x1, x2]                                
		3			(-infty, +infty)
		4		(-infty, x1] U [x2, infty)
		5		(-infty, x1] U [x2, x3] U [x4, infty)
		6		    [x1, x2] U [x3, x4]

-----------------------------------------------------------------------
NB: BPP modified the type codes used by AR and LM in the prequel to 
    match the ones used by CLR.
*/
		
		tempname cross zpz sqrtzpzi zpy MM ypz sqrtomegai v 	/*
			*/	d M N alpha C A D aa x1 x2 g type
		if ("`instcons'" == "") {
			mat accum `cross' = `inst' `one' `y1' `y2' 	/*
				*/	if `touse', noconstant
		}
		else {
			mat accum `cross' = `inst' `y1' `y2' 		/*
				*/	if `touse', noconstant
		}
	  
		mat `zpz' = `cross'[1..`k', 1..`k']
		mat `zpy' = `cross'[1..`k', (`k'+1)...]
		mat_inv_sqrt `zpz' `sqrtzpzi'
		mat_inv_sqrt `omega' `sqrtomegai'
		mat `ypz'=`zpy''
		mat `MM' = `sqrtomegai'*`ypz'*inv(`zpz')*`zpy'*`sqrtomegai'
		loc `k'=`df'
		mat symeigen `v' `d' = `MM'
		sca `M' = `d'[1,1]
		sca `N' =`d'[1,2]
		sca `alpha' = 1-`level'/100

/* inversion of CLR*/
		inversefun `M' `df' `alpha' `C'

		mat `A' =inv(`omega')*`ypz'*inv(`zpz')*`zpy'*inv(`omega')- /*
			*/	`C'*inv(`omega')
		sca `D' = -det(`A')
		sca `aa' = `A'[1,1]
		if (`aa'<0) {
	 		if (`D' <0) {
				sca `type'=1
			}
	 		else{
				sca `type'=2
				sca `x1'= (-`A'[1,2] + sqrt(`D'))/`aa'
				sca `x2' = (-`A'[1,2] - sqrt(`D'))/`aa'
				mat `g'=(`x1'\ `x2')
 			}
 		}
		else{
			if (`D'<0) {
		  		sca `type'=3
			}
	 		else {
		  		sca `type'=4
		  		sca `x1'= (-`A'[1,2]-sqrt(`D'))/`aa'
		  		sca `x2'= (-`A'[1,2]+sqrt(`D'))/`aa'
		  		mat `g'=(`x1' \ `x2')
 			}
	 	}
	 	ereturn local LR_type `"`=`type''"'
	 	if `type' == 2 | `type' == 4 {
	 		eret scalar LR_x1 = `x1'
	 		eret scalar LR_x2 = `x2'
	 	}

/* inversion of LM */

if "`LM'" != ""	{		/* BPP add 20051220 */
		tempname lmcv q1 q2 A1 A2 D1 D2 y1 y2 y3 y4 type1
		sca `lmcv' = invchi2tail(1, (1-`level'/100))
		if (`df'==1) {
			sca `q1' = `M'-`lmcv'
			mat `A1' =inv(`omega')*`ypz'*inv(`zpz')*`zpy'*	/*
				*/	inv(`omega')-`q1'*inv(`omega')
			sca `D1' = -4*det(`A1')
			sca `y1'= (-2*`A1'[1,2]+ sqrt(`D1'))/2/`A1'[1,1]
			sca `y2'= (-2*`A1'[1,2]-sqrt(`D1'))/2/`A1'[1,1]
			if (`A1'[1,1]>0) { 
				if (`D1'>0) {
					sca `type1'=4 
						/* two infinite intervals*/
					eret scalar LM_x1 = `y2'
					eret scalar LM_x2 = `y1'
				}
				else {
					sca `type1'=3
				}
			}
			else{
				if (`D1'>0) {
					sca `type1'=2 /*one interval */
					eret scalar LM_x1 = `y1'
					eret scalar LM_x2 = `y2'

				}
				else {
					sca `type1'=3
				}
			}
		}
		else {
			if ((`M' +`N' - `lmcv')^2-4*`M'*`N'<0) { 
		 		sca `type1' = 3
		 	}
			else {
			    sca `q1' = (`M'+ `N' - `lmcv' - 	/*
	 			*/ sqrt((`M'+`N' - `lmcv')^2 -  /*
	 			*/	4*`M'*`N'))/2
	 		    sca `q2' = (`M'+`N' - `lmcv' + 		/*
	 			*/ sqrt((`M'+`N'-`lmcv')^2 - 	/*
	 			*/	4*`M'*`N'))/2
	 		    if ((`q1' < `N') | (`q2' > `M')) {
	 			sca `type1' = 3
	 		    }
	 		    else { 		
				mat `A1' = inv(`omega')*`ypz'*inv(`zpz')*   /*
					*/ `zpy'*inv(`omega')-`q1'*inv(`omega')
				mat `A2' = inv(`omega')*`ypz'*inv(`zpz')*   /*
					*/ `zpy'*inv(`omega')-`q2'*inv(`omega')
		 		sca `D1' = -4*det(`A1')
				sca `D2' = -4*det(`A2')
	 			if (`A1'[1,1]>0) { 
		  			if (`A2'[1,1]>0) { 
						sca `type1' = 5
						sca `y1' = (-2*`A1'[1,2] +  /*
						    */ sqrt(`D1'))/2/`A1'[1,1]
						sca `y2' = (-2*`A1'[1,2] -  /*
						    */ sqrt(`D1'))/2/`A1'[1,1]
						sca `y3' = (-2*`A2'[1,2] +  /*
						    */ sqrt(`D2'))/2/`A2'[1,1]
						sca `y4' = (-2*`A2'[1,2] -  /*
						    */ sqrt(`D2'))/2/`A2'[1,1]
						eret scalar LM_x1 = `y4'
						eret scalar LM_x2 = `y2'
						eret scalar LM_x3 = `y1'
						eret scalar LM_x4 = `y3'
					}
		  			else {
						sca `type1' = 6
						sca `y1' = (-2*`A1'[1,2] +  /*
						    */ sqrt(`D1'))/2/`A1'[1,1]
						sca `y2' = (-2*`A1'[1,2] -  /*
						    */ sqrt(`D1'))/2/`A1'[1,1]
						sca `y3' = (-2*`A2'[1,2] +  /*
						    */ sqrt(`D2'))/2/`A2'[1,1]
						sca `y4' = (-2*`A2'[1,2] -  /*
						    */ sqrt(`D2'))/2/`A2'[1,1]
						eret scalar LM_x1 = `y3'
						eret scalar LM_x2 = `y4'
						eret scalar LM_x3 = `y2'
						eret scalar LM_x4 = `y1'
			  		}
			  	}
				if (`A1'[1,1]<=0) {
					sca `type1' =5
			  		sca `y1' = (-2*`A1'[1,2] + 	/*
			  			*/  sqrt(`D1'))/2/`A1'[1,1]
					sca `y2' = (-2*`A1'[1,2] - 	/*
						*/  sqrt(`D1'))/2/`A1'[1,1]
					sca `y3' = (-2*`A2'[1,2] + 	/*
						*/  sqrt(`D2'))/2/`A2'[1,1]
					sca `y4' = (-2*`A2'[1,2] -	/*
						*/  sqrt(`D2'))/2/`A2'[1,1]
					eret scalar LM_x1 = `y1'
					eret scalar LM_x2 = `y3'
					eret scalar LM_x3 = `y4'
					eret scalar LM_x4 = `y2'
		  		}
			    }
			}
		}
		eret local LM_type `"`=`type1''"'
}	/* End of LM block */

/* inversion of AR */
if "`AR'" != ""	{		/* BPP add 20051220 */
		tempname lmcv1  AAA type2 xx1 xx2 DDD aaa
		sca `lmcv1' = invchi2tail(`df', (1-`level'/100))
		mat `AAA' =`ypz'*inv(`zpz')*`zpy'-`lmcv1'*`omega'
		sca `DDD' = -det(`AAA')
		sca `aaa' = `AAA'[2,2]
		if (`aaa'<0) {
	 		if (`DDD' <0) {
				sca `type2'=3
			}
	 		else{
				sca `type2'=4
		 		sca `xx1'= (`AAA'[1,2] + sqrt(`DDD'))/`aaa'
		 		sca `xx2' = (`AAA'[1,2] - sqrt(`DDD'))/`aaa'
				eret scalar AR_x1 = `xx1'
				eret scalar AR_x2 = `xx2'
 			 }
		}
		else {
			if (`DDD'<0) {
				sca `type2'=1
			}
	 		else {
		  		sca `type2'=2
		  		sca `xx1'= (`AAA'[1,2]-sqrt(`DDD'))/`aaa'
		  		sca `xx2'= (`AAA'[1,2]+sqrt(`DDD'))/`aaa'
				eret scalar AR_x1 = `xx1'
				eret scalar AR_x2 = `xx2'
 			}
	 	}
		eret local AR_type `"`=`type2''"'
}	/* End of AR block */

	} /*---end of quiet--*/

/*------Testing block-----*/
	tempname  a b aprime bprime	
	
	/* Set up a and b vectors.	*/
	mat `a' = (`bbb'\1)
	mat `b' = (1\(-1*`bbb'))
	mat `aprime' = `a''		/* Macro parser chokes 	*/
	mat `bprime' = `b''		/* otherwise.		*/

 	tempname oia
	mat `oia' = inv(`omega')*`a'
	/* Compute scalars b'*O*b and a'*O^-1*a	*/
	tempname bob aoia
	matrix `bob' = `bprime'*`omega'*`b'
	scalar `bob' = trace(`bob')
	matrix `aoia' = `aprime'*inv(`omega')*`a'
	scalar `aoia' = trace(`aoia')
	
	tempname  sbar tbar
	mat `sbar' = `sqrtzpzi'*`zpy'*`b'/sqrt(`bob')
	mat `tbar' = `sqrtzpzi'*`zpy'*`oia'/sqrt(`aoia')
 
	tempname ar lm lr wald qt pval_new
	calcstat1 `sbar' `tbar' `ar' `lm' `qt'
	calcstat2 `sbar' `tbar' `aoia' `bob' `bbb' `omega' `lr' 
	loc kk=`k'-1

	if (("`instcons'" == "")&("`cons'" == "")) {
		new_try `kk' `qt' `lr' `pval_new'
	}
	else {
		new_try `k' `qt' `lr' `pval_new'
	}
		
	if "`AR'" != "" {
		tempname arcv arpv 
		if (("`instcons'" == "")&("`cons'" == "")) {
			sca `arcv' = invchi2tail((`k'-1), (1-`level'/100))
			sca `arpv' = 1-chi2((`k'-1), `ar')
		}
		else {
			sca `arcv' = invchi2tail(`k', (1-`level'/100))
			sca `arpv' = 1-chi2(`k', `ar')	
		}
		eret scalar p_AR = `arpv'
	}
	if "`LM'" != "" {
		tempname lmpv
		sca `lmpv' = 1-chi2(1, `lm')
		eret scalar p_LM = `lmpv'
	}
	
	eret scalar p_LR = `pval_new'

	eret local cmd "condivreg"

	myprint		/* Level stored in e(level) */
	myprintregions "`interval'" "`AR'" "`LM'"
		
end
  
	

/*----------------------------------------------------------------------*/

/* Computes AR and LM statistics.	*/
prog def calcstat1

	args sbar tbar ar lm qt
	tempname sbarp tbarp
	
	mat `sbarp' = `sbar''
	mat `tbarp' = `tbar''	
	mat `ar' = `sbarp'*`sbar'
	sca `ar' = trace(`ar')
	mat `lm' = (trace(`sbarp'*`tbar')^2) / trace(`tbarp'*`tbar')
	sca `lm' = trace(`lm')
	mat `qt' = `tbarp'*`tbar'
	sca `qt' = trace(`qt')

	
end 

/* Computes LR and Wald statistics.	*/
prog def calcstat2

	args sbar tbar aoia bob beta omega lr 
	tempname sbarp tbarp
	
	mat `sbarp' = `sbar''
	mat `tbarp' = `tbar''		
	tempname ss st tt
	mat `ss' = `sbarp'*`sbar'
	sca `ss' = trace(`ss')
	mat `tt' = `tbarp'*`tbar'
	sca `tt' = trace(`tt')
	mat `st' = `sbarp'*`tbar'
	sca `st' = trace(`st')
	sca `lr' = 0.5*(`ss' - `tt' + sqrt((`ss' + `tt')^2 -	/*
			*/	4*(`ss'*`tt' - (`st')^2)))	
  
end



program define new_try

	args k qt lrstat pval_new

	tempname gamma pval  u s2 qs farg1 farg2 farg wt
	
	sca `gamma' = 2*exp(lngamma(`k'/2)) / sqrt(_pi) / exp(lngamma((`k'-1)/2))

	if("`k'" == "1") {
		sca `pval' = 1 - chi2(`k', `lrstat')
	}
	else if ("`k'"== "2") {
		local ni 20
		mat `u' = J(`ni'+1,1,0)
		mat `s2' = J(`ni'+1,1,0)
		mat `qs' = J(`ni'+1,1,0)
		mat `wt' = J(1,`ni'+1,2)
		mat `farg1' = J(`ni'+1,1,0)
		mat `qs'[1,1] = (`qt'+`lrstat')
		mat `farg1'[1,1] = `gamma'*chi2(`k',`qs'[1,1])
		forv i =1(1)`ni'{
			mat `u'[`i'+1,1] = `i'*_pi/2/`ni'
			mat `s2'[`i'+1,1] = sin(`u'[`i'+1,1])
			mat `qs'[`i'+1,1] = (`qt'+`lrstat') / 	/*
				*/ (1+(`qt'/`lrstat')*`s2'[`i'+1,1]*`s2'[`i'+1,1])
			mat `farg1'[`i'+1,1] = `gamma'*chi2(`k',`qs'[`i'+1,1])
		}
		mat `wt'[1,1] = 1
		mat `wt'[1,`ni'+1] = 1
		local ni = `ni'/2
		forv i =1(1)`ni'{
			mat `wt'[1,`i'*2] = 4
		}
		local ni = `ni'*2
		mat `wt' = `wt'*_pi/2/3/`ni'
		mat `pval' = `wt'*`farg1'
		sca `pval' = 1-trace(`pval')
	}
	else if ("`k'"== "3") {
		local ni 20
		mat `s2' = J(`ni'+1,1,0)
		mat `qs' = J(`ni'+1,1,0)
		mat `wt' = J(1,`ni'+1,2)
		mat `farg1' = J(`ni'+1,1,0)
		mat `qs'[1,1] = (`qt'+`lrstat')
		mat `farg1'[1,1] = `gamma'*chi2(`k',`qs'[1,1])
		forv i =1(1)`ni'{
			mat `s2'[`i'+1,1] = `i'/`ni'
			mat `qs'[`i'+1,1] = (`qt'+`lrstat') / 	/*
				*/ (1+(`qt'/`lrstat')*`s2'[`i'+1,1]*`s2'[`i'+1,1])
			mat `farg1'[`i'+1,1] = `gamma'*chi2(`k',`qs'[`i'+1,1])
		}
		mat `wt'[1,1] = 1
		mat `wt'[1,`ni'+1] = 1
		local ni = `ni'/2
		forv i =1(1)`ni'{
			mat `wt'[1,`i'*2] = 4
		}
		local ni = `ni'*2
		mat `wt' = `wt'/3/`ni'
		mat `pval' = `wt'*`farg1'
		sca `pval' = 1-trace(`pval')
	}
	else if ("`k'"== "4") {
		local eps .02
		local ni 50
		mat `s2' = J(`ni'+1,1,0)
		mat `qs' = J(`ni'+1,1,0)
		mat `wt' = J(1,`ni'+1,2)
		mat `farg' = J(`ni'+1,1,0)
		mat `farg1' = J(`ni'+1,1,0)
		mat `farg2' = J(`ni'+1,1,1)
		mat `qs'[1,1] = (`qt'+`lrstat')
		mat `farg1'[1,1] = `gamma'*chi2(`k',`qs'[1,1])
		mat `farg'[1,1] = `farg1'[1,1]*`farg2'[1,1]
		forv i = 1(1)`ni'{
			mat `s2'[`i'+1,1] = `i'/`ni'*(1-`eps')
			mat `qs'[`i'+1,1] = (`qt'+`lrstat') / 	/*
				*/ (1+(`qt'/`lrstat')*`s2'[`i'+1,1]*`s2'[`i'+1,1])
			mat `farg1'[`i'+1,1] = `gamma'*chi2(`k',`qs'[`i'+1,1])
			mat `farg2'[`i'+1,1] = sqrt(1-`s2'[`i'+1,1]*`s2'[`i'+1,1])
			mat `farg'[`i'+1,1] = `farg1'[`i'+1,1]*`farg2'[`i'+1,1]
		}
		mat `wt'[1,1] = 1
		mat `wt'[1,`ni'+1] = 1
		local ni = `ni'/2
		forv i = 1(1)`ni'{
			mat `wt'[1,`i'*2] = 4
		}
		local ni = `ni'*2
		mat `wt' = `wt'/3/`ni'*(1-`eps')
		mat `pval' = `wt'*`farg'
		sca `pval' = 1-trace(`pval')
		sca `s2' = 1-`eps'/2
		sca `qs' = (`qt'+`lrstat')/(1+(`qt'/`lrstat')*`s2'*`s2')
		sca `farg1' = `gamma'*chi2(`k',`qs')
		sca `farg2' = 0.5*(asin(1)-asin(1-`eps'))-(1-`eps') / /*
			*/	2*sqrt(1-(1-`eps')*(1-`eps'))
		sca `pval' = `pval'-`farg1'*`farg2'
	}
	else {
		local ni 20
		mat `s2' = J(`ni'+1,1,0)
		mat `qs' = J(`ni'+1,1,0)
		mat `wt' = J(1,`ni'+1,2)
		mat `farg' = J(`ni'+1,1,0)
		mat `farg1' = J(`ni'+1,1,0)
		mat `farg2' = J(`ni'+1,1,1)
		mat `qs'[1,1] = (`qt'+`lrstat')
		mat `farg1'[1,1] = `gamma'*chi2(`k',`qs'[1,1])
		mat `farg'[1,1] = `farg1'[1,1]*`farg2'[1,1]
		forv i =1(1)`ni'{
			mat `s2'[`i'+1,1] = `i'/`ni'
			mat `qs'[`i'+1,1] = (`qt'+`lrstat') / 	/*
				*/ (1+(`qt'/`lrstat')*`s2'[`i'+1,1]*`s2'[`i'+1,1])
			mat `farg1'[`i'+1,1] = `gamma'*chi2(`k',`qs'[`i'+1,1])
			if "`i'" == "`ni'"{
				mat `farg2'[`i'+1,1] = 0
			}
			else{
				mat `farg2'[`i'+1,1] = /*
				*/   (1-`s2'[`i'+1,1]*`s2'[`i'+1,1])^((`k'-3)/2)
			}
			mat `farg'[`i'+1,1] = `farg1'[`i'+1,1]*`farg2'[`i'+1,1]
		}
		mat `wt'[1,1] = 1
		mat `wt'[1,`ni'+1] = 1
		local ni = `ni'/2
		forv i = 1(1)`ni'{
			mat `wt'[1,`i'*2] = 4
		}
		local ni = `ni'*2
		mat `wt' = `wt'/3/`ni'
		mat `pval' = `wt'*`farg'
		sca `pval' = 1-trace(`pval')
	}
	
	sca `pval_new' = `pval'

end 


prog def inversefun

	args M k alpha C

	tempname eps a  b x fa fb lrstat fx

	sca `eps' = 0.000001
	sca `a' = `eps' 
	sca `b' = `M' - `eps'
	sca `lrstat'= `M' - `a'

	new_try `k' `a' `lrstat' `fa'
	sca `lrstat' = `M' - `b'
	new_try `k' `b' `lrstat' `fb'

	if(`fa' > `alpha') {
		sca `C' = `a'
	}
	else  if ( `fb' <`alpha') {
		sca `C' = `b'
	}
	else {
		while (`b'-`a'>`eps') {
			sca `x' = (`b'-`a')/2+`a'
			sca `lrstat'= `M'-`x'
			new_try `k' `x' `lrstat' `fx'
			if (`fx' >`alpha') {
			 	  sca `b' = `x'
			}
			else {
				sca `a' = `x'
			}				 
		}
		sca `C' = `x'
	}
	  
end



prog def mat_inv_sqrt

	args in out
	tempname v vpri lam srlam

	loc k = rowsof(`in')
	mat symeigen `v' `lam' = `in'
	mat `vpri' = `v''
	/* Get sqrt(lam)	  */
	mat `srlam' = diag(`lam')
	forv i = 1/`k' {
		mat `srlam'[`i', `i'] = 1/sqrt(`srlam'[`i', `i'])
	}
	mat `out' = `v'*`srlam'*`vpri'

end

program define my2sls, rclass

	args ry1 ry2 rexog rinst cons instcons touse

	quietly {
		tempname b2sls var2sls rss rssold
		tempvar y2hat
		if "`instcons'" == "" {
			reg `ry2' `rexog' `rinst' if `touse'
		}
		else {
			reg `ry2' `rexog' `rinst' if `touse', noconstant
		}
		predict double `y2hat' if `touse', xb
		reg `ry1' `y2hat' `rexog' if `touse', `cons'
		mat `b2sls' = e(b)
		mat `var2sls' = e(V)
		sca `rssold' = e(rss)
		tempvar resids hat
		replace `y2hat' = `ry2'
		predict double `resids' if `touse', residuals
		replace `resids' = `resids'^2
		su `resids', meanonly
		sca `rss' = r(sum)
		mat `var2sls' = `var2sls'*`rss'/`rssold'
		return scalar rss = `rss'
		return matrix beta `b2sls'
		return matrix var `var2sls'
	}

end

program define myliml, rclass

	args ry1 ry2 rexog rinst cons instcons touse

	quietly {
		/* LIML estimator.	*/
		/* Follows Davidson and MacKinnon (1993).	*/
		tempvar one
		gen double `one' = 1
		if "`instcons'" == "" {
			loc x2 `rinst' `one'
		}
		else {
			loc x2 `rinst'
		}
		if "`cons'" == "" {
			loc x1 `rexog' `one'
		}
		else {
			loc x1 `rexog'
		}
	
		tempname xpx x1px1 x1py1 y1py1
		mat accum `xpx' = `ry2' `x1' if `touse', noconstant
		mat `y1py1' = `xpx'[1, 1]

		/* Compute Y1'MxY1.	*/
		tempname y1mxy1
		loc k4 : word count "`ry2'"
		if `k4' != 1 {
			di as error /*
			*/ "Multiple endogenous RHS variables not implemented."
			exit 198
		}
		reg `ry2' `x1' `x2' if `touse', noconstant 
			/* May be two ones here, but reg will catch it. */
		tempvar y2resid
		predict double `y2resid' if `touse', residuals
		replace `y2resid' = `y2resid'^2
		summ `y2resid', meanonly
		sca `y1mxy1' = r(sum)
	
		/* Compute kappa.	*/
		/* First get Y'MxY.	*/
		tempvar y1h y2h y1hy2h
		tempname s11 s12 s22 ymxy ym1y
		reg `ry1' `x1' `x2' if `touse', noconstant
		predict double `y1h' if `touse', residuals
		reg `ry2' `x1' `x2' if `touse', noconstant
		predict double `y2h' if `touse', residuals
		gen double `y1hy2h' = `y1h'*`y2h' if `touse'
		replace `y1h' = `y1h'^2
		replace `y2h' = `y2h'^2
		summ `y1hy2h' if `touse', meanonly
		scalar `s12' = r(sum)
		summ `y1h' if `touse', meanonly
		scalar `s11' = r(sum)
		summ `y2h' if `touse', meanonly
		scalar `s22' = r(sum)
		mat `ymxy' = J(2,2,0)
		mat `ymxy'[1, 1] = `s11'
		mat `ymxy'[1, 2] = `s12'
		mat `ymxy'[2, 1] = `s12'
		mat `ymxy'[2, 2] = `s22'

		/* And Y'M1Y.	*/
		tempvar y1hh y2hh y1hhy2hh
		reg `ry1' `x1' if `touse', noconstant
		predict double `y1hh' if `touse', residuals
		reg `ry2' `x1' if `touse', noconstant
		predict double `y2hh' if `touse', residuals
		gen double `y1hhy2hh' = `y1hh'*`y2hh' if `touse'
		replace `y1hh' = `y1hh'^2
		replace `y2hh' = `y2hh'^2
		summ `y1hhy2hh' if `touse', meanonly
		scalar `s12' = r(sum)
		summ `y1hh' if `touse', meanonly
		scalar `s11' = r(sum)
		summ `y2hh' if `touse', meanonly
		scalar `s22' = r(sum)
		mat `ym1y' = J(2,2,0)
		mat `ym1y'[1, 1] = `s11'
		mat `ym1y'[1, 2] = `s12'
		mat `ym1y'[2, 1] = `s12'
		mat `ym1y'[2, 2] = `s22'

		/* Now form the matrix to get kappa from. */
		tempname ymxyih kmat vals vecs kappa bliml varliml
		mat_inv_sqrt `ymxy' `ymxyih'
		mat `kmat' = `ymxyih'*`ym1y'*`ymxyih'
		mat symeigen `vecs' `vals' = `kmat'
		loc junk = colsof(`vals')
		scalar `kappa' = `vals'[1, `junk'] 
		
		/* Now assemble the "X'X" matrix.	*/
		tempname xpxi xpy xpy1 xpy2 bliml varliml rss sighatsq
		mat `xpx'[1, 1] = `y1py1' - `kappa'*`y1mxy1'
		mat `xpxi' = syminv(`xpx')
		/* And the "X'y" matrix.	*/
		mat accum `xpy1' = `ry2' `ry1' if `touse', noconstant
		mat `xpy1' = `xpy1'[1, 2]
		mat `xpy1' = `xpy1' - `kappa'*`ymxy'[1, 2]
		mat accum `xpy2' = `x1' `ry1' if `touse', noconstant
		loc k1 : word count `x1'
		loc k2 = `k1' + 1
		mat `xpy2' = `xpy2'[1..`k1', `k2'...]
		mat `xpy' = `xpy1' \ `xpy2'
		mat `bliml' = (`xpxi'*`xpy')'
		
		/* Now compute sigma_squared to get the covariance matrix. */
		tempvar linpred residsq
		mat score `linpred' = `bliml' if `touse'
		gen double `residsq' = (`ry1' - `linpred')^2 if `touse'
		su `residsq' if `touse', meanonly
		sca `rss' = r(sum)
		sca `sighatsq' = r(sum) / r(N)
		mat `varliml' = `xpxi' * `sighatsq'

		return matrix beta `bliml'
		return matrix var `varliml'
		return scalar rss = `rss'
	}	
	
end	

prog def myprint

	local level `e(level)'
	di
	if "`e(model)'" == "2sls" {
		di as text "Instrumental variables (2SLS) regression"
	}
	else {
		di as text "Instrumental variables (LIML) regression"
	}
	di
	di as text "First-stage results" _col(56) "Number of obs =" 	/*
		*/ as result %8.0f `e(N)'
	if "`e(model)'" == "2sls" {
		di as text "{hline 23}" _col(56) "F(" %3.0f `e(df_m)' 	/*
			*/ "," %6.0f `e(df_r)' ") =" as result %8.2f `e(F)'
		di as text "F(" %3.0f `e(df_m_first)' "," %6.0f 	/*
			*/ `e(df_r_first)' ") =" as result %8.2f 	/*
			*/ `e(F_first)' _col(56) as text "Prob > F      =" /*
			*/ as result %8.4f Ftail(`e(df_m)', `e(df_r)', `e(F)')
	}
	else {
		di as text "{hline 23}" _col(56) "Wald chi2(" %2.0f 	/*
			*/ `e(df_m)' ") =" as result %8.2f `e(wald)'
		di as text "F(" %3.0f `e(df_m_first)' "," %6.0f 	/*
			*/ `e(df_r_first)' ") =" as result %8.2f 	/*
			*/ `e(F_first)' _col(56) as text 		/*
			*/ "Prob > w      =" as result %8.4f 	/*
			*/ chi2tail(`e(df_m)', `e(wald)')
	}	
	di as text "Prob > F      =" as result %8.4f 		/*
		*/ Ftail(`e(df_m_first)', `e(df_r_first)', `e(F_first)') /*
		*/ _col(56) as text "R-squared     =" as result %8.4f `e(r2)'
	di as text "R-squared	  =" as result %8.4f `e(r2_first)'	/*
		*/ _col(56) as text "Adj R-squared =" as result %8.4f `e(r2_a)'
	di as text "Adj R-squared =" as result %8.4f `e(r2_a_first)'	/*
		*/ _col(56) as text "Root MSE      =" as result 	/*
		*/ %8.3f `e(rmse)' 
	di
	eret display, level(`level')
	di as text "Instrumented:  " _c
	Disp `e(instd)'
	di as text "Instruments:   " _c
	if "`e(instcons)'" == "yes" {
		Disp `e(insts)'
	}
	else {
		Disp `e(insts)' (No constant included)
	}
	di as text "Confidence set and p-value for " 			/*
		*/ abbrev("`e(instd)'", 13) " "				/*
		*/ "are based on normal approximation"
	di as text "{hline 78}"
	
	
end

program define myprintregions

	args interval AR LM

	di _n _n
	di as text "{hline 78}"
	if ("`interval'"=="") {
		local what "confidence set"
	}
	else {
		local what "confidence interval"
	}
	local plural
	if "`AR'`LM'" != "" {
		local plural "s"
	}
	local title : di "Coverage-corrected `what'`plural' and p-value`plural'"
	local titcol `=int((78 - length("`title'"))/2 + 1)'
	di as text _col(`titcol') "`title'"
	
	local title : di "for Ho: _b[`e(instd)'] = " %-9.0g e(H0_b)
	local titcol `=int((78 - length("`title'"))/2 + 1)'
	di as text _col(`titcol') "`title'"
	
	local title : di "LIML estimate of _b[`e(instd)'] = " /*
		*/	%-9.0g `e(b_instd_liml)'
	local titcol `=int((78 - length("`title'"))/2 + 1)'
	di as text _col(`titcol') "`title'"
	
	di
	di as text "{hline 78}"
	local title : di proper("`what'")
	local titcol `=int((70 - 15 - length("`title'"))/2) + 15'
	di as text _col(2) "Test" _col(`titcol') "`title'" _col(71) "p-value"
	di as text "{hline 78}"	
	
	local lrpval  "as result _col(72) %6.4f e(p_LR)"
	if `e(LR_type)' == 1 {
		local lrint "empty"
	}
	else if `e(LR_type)' == 2 {
		local lrint : di "[" 					/*
			*/	 %9.0g e(LR_x1)				/*
			*/	 ", "					/*
			*/	 %9.0g e(LR_x2)				/*
			*/	 "]"
	}
	else if `e(LR_type)' == 3 {
		local lrint "(-inf, +inf)"
	}
	else {
		if "`interval'" == "" {
			local lrint : di "(-inf, " 			/*
				*/	 %9.0g e(LR_x1) 		/*
				*/       "] U [" 			/*
				*/	 %9.0g e(LR_x2) 		/*
				*/	 ", +inf)" 
		}
		else {
			local lrint "(-inf, +inf)"
		}
	}
	local len `=length(`"`lrint'"')'
	local len `=int((70 - 15 - `len')/2)'
	local len `=`len'+15'
	di _col(2) "Conditional LR" _col(`len') as result "`lrint'" `lrpval'

	if ("`AR'"!="") {
		local arpval "as result _col(72) %6.4f e(p_AR)"
		if `e(AR_type)' == 1 {
			local arint "empty"
		}
		else if `e(AR_type)' == 2 {
			local arint : di "[" 				/*
				*/	 %9.0g e(AR_x1)			/*
				*/	 ", "				/*
				*/	 %9.0g e(AR_x2)			/*
				*/	 "]"
		}
		else if `e(AR_type)' == 3 {
			local arint "(-inf, +inf)"
		}
		else {
			if "`interval'" == "" {
				local arint : di "(-inf, " 		/*
					*/	 %9.0g e(AR_x1) 	/*
					*/       "] U [" 		/*
					*/	 %9.0g e(AR_x2) 	/*
					*/	 ", +inf)" 
			}
			else {
				local arint "(-inf, +inf)"
			}
		}
		local len `=length(`"`arint'"')'
		local len `=int((70 - 15 - `len')/2)'
		local len `=`len'+15'
		di _col(2) as text "Anderson-Rubin" 			/*
			*/	_col(`len') as result "`arint'" `arpval'
	}
	if ("`LM'"!="") {
		local lmpval "as result _col(72) %6.4f e(p_LM)"
		if `e(LM_type)' == 2 {
			local lmint : di "["				/*
				*/	 %9.0g e(LM_x1)			/*
				*/	 ", "				/*
				*/	 %9.0g e(LM_x2)			/*
				*/	"]"
		}
		else if `e(LM_type)' == 3 {
			local lmint "(-inf, +inf)"
		}
		else if `e(LM_type)' == 4 {
			if "`interval'" == "" {
				local lmint : di "(-inf, " 		/*
					*/	 %9.0g e(LM_x1) 	/*
					*/       "] U [" 		/*
					*/	 %9.0g e(LM_x2) 	/*
					*/	 ", +inf)" 
			}
			else {
				local lmint "(-inf, +inf)"
			}
		}
		else if `e(LM_type)' == 5 {
			if "`interval'" == "" {
				local lmint : di "(-inf, "		/*
					*/	 %6.0g e(LM_x1)		/*
					*/	 "] U ["		/*
					*/	 %6.0g e(LM_x2)		/*
					*/	 ", "			/*
					*/	 %6.0g e(LM_x3)		/*
					*/	 "] U ["		/*
					*/	 %6.0g e(LM_x4)		/*
					*/	 ", +inf)"	
			}
			else {
				local lmint "(-inf, +inf)"
			}
		}
		else if `e(LM_type)' == 6 & e(LM_x1) < e(LM_x3) {
			if "`interval'" == "" {
				local lmint : di "["			/*
					*/	 %9.0g e(LM_x1)		/*
					*/	 ", "			/*
					*/	 %9.0g e(LM_x2)		/*
					*/	 "] U ["		/*
					*/	 %9.0g e(LM_x3)		/*
					*/	 ", "			/*
					*/	 %9.0g e(LM_x4)		/*
					*/	 "]"
			}
			else {
				local lmint : di "["			/*
					*/	 %9.0g e(LM_x1)		/*
					*/	 ", "			/*
					*/	 %9.0g e(LM_x4)		/*
					*/	 "]"
			}
		}
		else {
			if "`interval'" == "" {
				local lmint : di "["			/*
					*/	 %9.0g e(LM_x3)		/*
					*/	 ", "			/*
					*/	 %9.0g e(LM_x4)		/*
					*/	 "] U ["		/*
					*/	 %9.0g e(LM_x1)		/*
					*/	 ", "			/*
					*/	 %9.0g e(LM_x2)		/*
					*/	 "]"			
			}
			else {
				local lmint : di "["			/*
					*/	 %9.0g e(LM_x3)		/*
					*/	 ", "			/*
					*/	 %9.0g e(LM_x2)		/*
					*/	 "]"
			}
		}
		local len `=length(`"`lmint'"')'
		local len `=int((70 - 15 - `len')/2)'
		local len `=`len'+15'
		di _col(2) as text "Score (LM)" 			/*
			*/ _col(`len') as result "`lmint'" `lmpval'
	}
	di as text "{hline 78}"

end


/* Lifted straight out of ivreg.ado.	*/
program define IsStop, sclass

	/* sic, must do tests one-at-a-time, 0 may be very large */
	if `"`0'"' == "[" {
		sret local stop 1
		exit
	}
	if `"`0'"' == "," {
		sret local stop 1
		exit
	}
	if `"`0'"' == "if" {
		sret local stop 1
		exit
	}
	if `"`0'"' == "in" {
		sret local stop 1
		exit
	}
	if `"`0'"' == "" {
		 sret local stop 1
		 exit
	}
	else	 sret local stop 0

end

/* More code borrowed from ivreg.ado.	*/
program define Disp
	local first ""
	local piece : piece 1 64 of `"`0'"'
	local i 1
	while "`piece'" != "" {
		di in gr "`first'`piece'"
		local first "              "
		local i = `i' + 1
		local piece : piece `i' 64 of `"`0'"'
	}
	if `i'==1 { 
		di 
	}
end

exit

